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Abstract 



In this paper we present a preliminary report on our work on the tracking of inter- 
nal layers in a singularly-perturbed convection-diffusion equation. We show why such 
tracking may be desirable, and we also show how to do it using domain decomposition 
based on asymptotic analysis. 


1 Introduction 


In this paper we present the analogue of a shock-tracking scheme for the resolution of an 
internal layer and its interaction with an ordinary boundary layer at the outflow. In the 
computation of compressible flows at high Mach number there has long been competition 
between shock tracking and shock capturing, and it is now generally agreed that both are 
needed. We generally And that the number of strong shocks is small, and they should be 
tracked in order to assure accuracy of the solution. For reasons of efficiency, however, the 
large number of weak shocks reverberating around the domain should be computed by a 
reliable shock-capturing scheme such as Roe’s method [14] or the method of Colella and 
Woodward [7]. Shocks are not always the most important features in fluid flows, but the 
tracking of such other phenomena is still far behind shock tracking. We first show why it 
may be desirable to track an internal layer, and then we show how such tracking may be 
accomplished via domain decomposition. 

For the sake of simplicity we consider a specific time- in dependent, singularly-perturbed, 
convection- diffusion equation 


a d x u -f- b d y u — e A u -f F 


( 1 . 1 ) 
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on a bounded domain ft in the (x,y)-plane. Here, A denotes the Laplace operator, and e 
is a small, positive number. For the moment, we impose Dirichlet conditions u = / on the 
boundary <9ft, but later we sometimes treat mixed boundary conditions. The function / 
is required to be piecewise smooth. (We use the term ‘smooth’ to mean some convenient 
degree of differentiability, say C 2 .) We assume that the coefficients a and 6 are smooth 
functions of x and y on fl. The source term F is assumed to be a smooth function of 
x } y, and u. Furthermore, we impose the restriction that |a| -f- |6| ^ 0 in the closure of 
0. This assumption implies that there are no stagnation points, and it greatly diminishes 
the complexity of the domain decomposition. Our assumption of semilinearity is much less 
restrictive because nonlinear problems are often solved via a sequence of linear problems 
with variable coefficients. Our discussion does not pertain to shock layers, however, since 
they violate the assumption of smoothness of a and b. , 

Previous work [4], [11], [12] on algorithms for (1.1) using domain decomposition based 
on asymptotic analysis has treated the special case of 6 = 0. It is true that a transformation 
of coordinates may be used to convert (1.1) to the case 6 = 0. In this paper we show why 
such a transformation is very desirable, and we present an algorithm to carry it out. 

The development of numerical methods for (1.1) in the singularly-perturbed case requires 
an understanding of the asymptotic behavior of its solutions as e J 0. We therefore begin 
with a brief description of the relationship between u and the solution U of the reduced 
equation 

ad x U + bdyU = F. ( 1 . 2 ) 

For more details see the books by Chang and Howes [1] and Eckhaus [8]. Equation (1.2) is 
easily solved by the method of characteristics, 


dx 





(1.3) 


The first two equations in (1.3) define characteristic curves. It is clear that we cannot 
impose the boundary condition U — f at every intersection of a characteristic curve with dfl. 
Instead, we subdivide dft into three sets, depending on the direction of the vector (a, 6). The 
inflow boundary Tj is the subset of dCl on which (a, 6) points into ft, the outflow boundary T 0 
is the subset of dfl on which (a, 6) points out of ft, and the tangential boundary Tt is the 
subset of dfl on which (a, 6) is tangent to dfl. For (1.3) the boundary condition U = / is 
imposed only on the inflow boundary T/. 


It is reasonable to expect to have u U for the solutions u of (1.1) wherever A u is not 
too large, i.e, wherever u is smooth. Because of the smoothness of the source term F and 
of the coefficients a and 6, the only mechanism for introducing nonsmooth behavior into the 
solution u of (1.1) is through the boundary condition u — f. One obvious difficulty is that 
we cannot force U = / on the outflow and tangential boundaries To U IV . The resolution of 
this difficulty is that there are boundary layers across which u changes rapidly from u K U 
to u = /. More precisely, when / is smooth the relation u & U is true except in the following 
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portions of SI. There may be what are called parabolic boundary layers along the tangential 
boundary IY, and there may be ordinary boundary layers along the outflow boundary Tq. 

Let us take a moment to explain the terminology ‘ordinary boundary layer’ and to de- 
scribe its properties. Consider a point P on To- In the vicinity of P we may construct a 
transformation (tr, r) > (a :,y) such that the origin (<7,t) = (0,0) is mapped to the point P. 
We may further require that the portion of a neighborhood of the origin with a > 0 is mapped 
into fl, so that a segment of the axis a = 0 is mapped onto a segment of the boundary Tq. 
In terms of the variables a and r the differential equation (1.1) takes the form 

« 

a d„u + bd T u = e(ci dlu + c 3 d a d T u + c 3 d\u) + F. (1.4) 

Note that the definition of outflow boundary implies that if c 4 is chosen to be positive, then 
it follows that a < 0. We have seen that we expect to have u & U away from the boundary 
<7 — 0 t while we require that u = f on the boundary <r — 0. That is, we expect u to vary 
slowly with respect to r but rapidly with respect to <r. Let us therefore introduce the scaling 

<j = eir, r = f (1.5) 

into (1.4). If we formally discard all but terms of the order of 1/e, we obtain a reduced 
equation 

ad s V = c x dlV. (1.6) 

The term ordinary boundary layer derives from the fact that (1.6) is an ordinary differential 
equation. The term exponential boundary layer is also used, because the solution V of (1.6) 
is the sum of a constant and an exponential function. Note that because 2/ci < 0, this 
exponential decreases with increasing a. Note also that in terms of the variable <r the rate of 
decrease is of the order of exp {— k<t/c}, where n is an average value of |a|/ci. We therefore 
expect the width of the ordinary boundary layer to be 0(c). . The book by Chang and 
Howes [1] gives theoretical justification for all of these heuristic manipulations. 

In the vicinity of the tangential boundary IY we use the characteristic curves (1.3) to 
define one set of coordinate lines, and we use them as the foundation for a local mapping 
{s,t) i — ► (a;,J/) in the vicinity of a point P on IY. In terms of these coordinates (1.1) takes 
the form 

d,u = c(c 4 d]u + c B d s d t u + c 6 d^u) + F. (1.7) 

We remark that the definition of flow direction implies that c$ > 0. We may require that a 
segment of the axis t = 0 maps onto a segment of the boundary IY and that positive values 
of t correspond to points in the interior of SI. Thus, the boundary layer has to accommodate 
a rapid transition from u « U for t > 0 to u ■ == / at t = 0. Let us therefore introduce the 
scaling 

3 = 5, t = T-Je (1.8) 

into (1.7). If we formally delete all terms smaller than 0(1), we obtain the reduced equation 

d i W = c 4 d?W + F. (1.9) 
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This equation is parabolic, giving rise to the term parabolic boundary layer. Furthermore, 
the thickness of the boundary layer for (1.9) is 0(1) with respect to t. With respect to t 
the parabolic boundary layer is therefore of width 0(y/i). Again, theoretical justification for 
these remarks may be found in [1]. 

If / has a discontinuity at a point P on Tj, then by (1.3) the function U will have a 
corresponding discontinuity in ft along the characteristic curve 7 through P. Similarly, if 
the Lie derivative of / along T/ is discontinuous at P, then grad U is discontinuous along 7. 
Because u is smooth, the lack of smoothness of U causes u to deviate substantially from U in 
a neighborhood of 7. As in the case of a parabolic boundary layer, if we introduce coordinates 
(s,<) derived from the characteristic variable s as given by (1.3), we find that (1.1) maps to 
an equation of the form (1.7) and that (1.9) is the appropriate reduced equation. We are 
therefore led to the conclusion that such an internal layer is parabolic in nature and that its 
width is 0(y/e). We again refer the reader to [1] for further justification. 

In the next section we analyze the behavior of a standard central difference scheme when 
there is an internal layer tilted at an angle to the grid, and we show that the numeri- 
cal approximation introduces downstream oscillations unless the internal layer is resolved. 
Therefore, we must use either a fine grid, an artificial increase of the viscosity e, or a grid 
aligned with the layer. Here we are discussing a grid effect, in [13] Hedstrom and Osterheld 
showed that the numerical errors for a coarse grid aligned with an internal layer are minimal 
even at large cell Reynolds numbers. 

In Section 3 we present an algorithm for the construction of an orthogonal grid with 
one coordinate direction aligned to the vector field (a, 6). This mapping requires the solu- 
tion of the telegraphers’ equation. In Section 4 we introduce a domain decomposition for 
a problem (1.1) with an internal layer and an ordinary boundary layer. In this domain de- 
composition the ordinary boundary layer and the internal layer have their own subdomains, 
and there is a separate subdomain for the region where these layers interact. In addition, 
in each subdomain a separate numerical method is used, depending on the local asymptotic 
behavior of the solution. 


2 A layer tilted to the grid. 

In this section we use a heuristic argument based on the modified equation to show why it is 
generally unwise to permit an internal layer not to be aligned with the grid. Specifically, we 
show that the standard central-difference scheme has grid effects which are modelled by a 
modified equation in the style of Warming and Hyett [16]. See Griffiths and Sanz-Serna [10] 
for a more modern exposition on modified equations. We shall see that the solutions of the 
modified equation are integrals of Airy functions, multiplied by a decaying exponential. The 
oscillations of this Airy function may or may not be completely damped by the exponential, 
depending on the values of a dimensionless parameter. We also derive the modified equation 
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for the upwind difference scheme, and as may be expected, we find that upwindhig adds 
numerical diffusion. 

For the discussion here we restrict our attention to the special case when the coefficients 
a and b in (1.1) are constant and the source term F vanishes. Then for convection with 
velocity V in the direction (cos 9, sin 0) we have the convection-diffusion equation 


V cos 9 d x u + V sin 9 d y u = e A u. (2.1) 

The reduced equation for (2.1) is 

V cos 9 d x U + V sin 9 d v U = 0, (2-2) 

and its characteristic curves are given by 

~ = V cos 9, ^ = V sin 9. (2.3) 

as is ' 

For the discussion here it suffices to restrict our attention to directions 0 < 9 < 7r/4. We 
remark that the special case 9 = 0 of flow parallel to the grid was examined by Hedstrom 
and Osterheld [13]. 

For (2.1) we use a rectangular domain 


0 = {(z>S/):0 < x < 1,0 < y < 1}. (2.4) 

Thus, under the conditions that 0 < 9 < 7r/4 the inflow boundary is at x — 0 and at y = 0, 
and the other two sides of the rectangle comprise the outflow boundary. On the inflow 
boundary we select a point of discontinuity y 0 > and we impose the conditions 


f 0.5( 1 + sgn(y - y 0 )) for x = 0, 
\ 0 for y = 0. 


(2.5) 


The discontinuity in the boundary data at yo induces an internal layer along the characteristic 
curve xsin# — (y — yo)cos9 = 0. In fact, the solution U to the reduced equation (2.2) with 
boundary data (2.5) is given by 


U = i i sgn(a: sin 9 - (y - y 0 ) cos 9). . (2.6) 

In order to minimize ordinary boundary layers along the outflow boundary, we impose the 
reduced equation (2.2) as boundary condition at x — 1 and at y — 1. * 


Consider the standard central-difference scheme for (2.1). We impose a square grid on Q 
with mesh size h, and we define the shift operators 


T x v(x,y) = v(x + h,y), T y v(x, y) = v(x,y + h). (2.7) 
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With this notation the central-difference approximation to (2.1) is 


V cos 9 

2h 


(T. 


T:')v + - T-')v = eDv, 


2h 


( 2 . 8 ) 


where D denotes the discrete Laplacian 


1 


D = j- 2 {T x + T y + T" 1 + T- 1 - 4/). 

On the inflow boundary T/ the boundary conditions for (2.8) are (2.5). On the outflow 
boundary To we use an upwind discretization of (2.2). 

* 

The modified equation for (2.8) is best written in terms of the rotated coordinate system 
aligned with the flow direction « 


s = x cos 9 + (y — t/ 0 ) sin 9 , 
t = —x sin 9 -f (y — y 0 ) cos 9. 


(2.9) 


We also introduce scalings of s and t in order to derive a modified equation in dimensionless 
form and to identify the pertinent parameters. It happens that for (2.1) or (2.8) on the 
halfplane x > 0 there is no natural length scale in the direction of the flow (the ^-direction). 
One may as well use a length scale L — 1. For the rectangular domain ft defined in (2.4) it 
is reasonable to take L to be the width of ft (L ~ 1) or the width of ft in the s-direction 
( L = sec#). We shall see that the natural scalings for the modified equation corresponding 
to the central difference method (2.8) are 


( 2 . 10 ) 


s — La, 

/Le \ ,/2 

t=T (v) ' 

Furthermore, the important dimensionless parameters for (2.8) are the cell Reynolds number 

Rh = ^- ( 2 . 11 ) 

and the degree of streamwise resolution 

h 

( 2 . 12 ) 

In terms of these parameters the modified equation for (2.8) is given by the following theorem. 


Theorem. Suppose that 0 < 9 < 7r/4. Suppose also that 7 <C 1 and that 7 <C Rh "C 
I/7. Then the modified equation for (2.8) is 


d'V = 8\v - + (j- + ^(3 + «*(«))) (2.13) 
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Remarks. The restriction that 7 = hj L <C 1 is reasonable for numerical computations, 
since we would want features in the streamwise direction to be resolved. The condition that 
7 <C Rh <C I/7 is also ordinarily satisfied in computations. We have written the modified 
equation in the form (2.13) in order to provide uniformity as sin(40) — » 0. The grid-induced 
oscillations appear only when sin(40) 7^ 0 and when Rh is moderately large. (Remember 
that we require iZ/,7 <C 1.) Under the condition that /Z^sin(4^) is bounded away from zero, 
the modified equation (2.13) takes the simpler form 


d ffV = d 2 T v - !^1 7 V2 jR 3/2 0 t V (2.14) 

In (2.14) the parameter 

r (2.i 5 ) 

measures the importance of the grid-induced numerical dispersion relative to the physical 
diffusion, and no numerically induced oscillations will be observed if /? is sufficiently small. 
For boundary data v = sgn(r) at r = 0 the solution of (2.14) may be expressed in terms of 
the Airy function, as is shown by Chin and Hedstrom [3]. In fact, a Fourier transformation 
with respect to r shows that 


v(<7, r) = — + f — — exp \if3cnJ 3 — <rw 2 + tVu>} doj. 

2 J — 00 llTTLO ^ ** 


(2.16) 


The reference [3] also provides figures and tables of the integrals (2.16) for various values 
of (3. The upshot is that whether or not there are oscillations depends on a parameter 


2a 1 / 3 

( 3 /?) 2 / 3 - 


(2.17) 


If a > 2, the diffusion is dominant, and there are no oscillations. For a < 2, however, there 
is a sequence of damped oscillations below the layer (r < 0). Because a is an increasing 
function of < 7 , as we proceed downwind a increases and the diffusion eventually removes the 
oscillations. 


Note with regard to the applicability of the modified equation that the internal layer is 
many grid cells wide and that the oscillations have wavelengths spanning many grid cells. 
This behavior makes the modified equation applicable, in that the derivation of a modified 
equation is based on the assumption that the numerical solution is smooth relative to the 
grid. The user of modified equations must always keep in mind that they are utterly useless 
in predicting variations in the numerical solution on the scale of 2 or 3 grid sizes or shorter. 

Finally, we also remark that the upwind difference scheme 

V — (/ - T~ l )u + V — (/ - T-> = c Du 

with the scalings (2.10) has the modified equation 

d a v = aj d 2 T v + a 2 d\v 
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with 


<*i = 1 + ^2 sin ( 2 ^) cos (f - o) , 

&2 = 4 (3 4 cos(40) 24(cos 3 9 4 sin 3 $)). 


The most significant numerical viscosity added by the grid effect is the deviation of c*i 
from 1. • 

Proof of the theorem. The idea of the modified equation is to make an ansatz that 
the solution of the difference scheme is smooth enough to be represented by a small number 
of terms of its Taylor expansion and to use this expansion to identify a partial differential 
equation which approximates (2.8) more closely than the original equation (2.1) does. 

Thus, we begin with the assumption that some smooth function v is a solution of (2.8). 
In this case the word ‘smooth’ is taken to mean that we may use Taylor approximations such 
as 

T x v=v + hd x v + ^ dlv 4y9’n| dfv (2.18) 

for the terms in (2.8). 

That is, we choose to neglect terms in the Taylor series approximation to (2.8) of order 
h 6 and higher. We therefore obtain the equation 


V cos 9 4 y- c? 3 e) 4 V sin 0 (d y v 4 ^ d 3 v) 

= e [d 2 x v 4 dlv 4 ^{dfv 4 dfv)) . 


The rotation of coordinates (2.9) then gives 

Vh 2 eh 2 

Vd.v 4 — L 3 [v] = e{d]v 4 dlv) 4 [*] 

with 


(2.19) 


L$\v\ — ^(3 4 cos(40)) dfv — | sin(40) d]d t v 4 | sin 2 (20) d,dfv 4 j sin(40) dfv , 

T 4 [u] = ^(3 4 cos(40))(c^u 4 dfv) — sin(40)(d 3 d t v — d,dfv) 4 3sin 2 (2 9)dfdfv. 
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Because we are interested in the effects of the internal layer, we expect derivatives of v 
in the i-direction (perpendicular to the layer) to be significantly larger than derivatives in 
the s-direction. The scalings (2.10) are selected to balance the terms d s v and dfv in (2.19). 
Thus, upon substituting (2.10) into (2.19) and dropping terms smaller than 0(^Rh)^ we 
obtain 

d 9 v = d 2 r v-(3 dlv + 1-dlv + 2^(3 + cos(40)) (2.20) 

Kh 4o 

with /? as defined by (2.15). 

The inexperienced user of modified equations may expect (2.20) to serve as a modified 
equation for (2.8). We cannot use it because the term involving d*v renders (2.20) unstable to 
high-frequency perturbations. The use of such a modified equation would predict numerical 
instabilities where there are none, and it is an instance of a modified equation not conforming 
to the difference scheme for phenomena of short wave length. The term d*v appears in (2.20) 
because we stopped the Taylor expansion (2.18) at d*v, and we went that far because the 
coefficient (3 of d*v in (2.20) is zero when sin(40) = 0. That is, we must replace d*v by 
something reasonable but harmless when /? is near zero, and for other values of (3 it need 
only be something harmless. Because d„v & d 2 v when /3 « 0, the substitution we make to 
render d*v harmless is that d*v as d 2 v. In this way a high-frequency instability is converted 
into an increase of dissipation in the streamwise direction, and it produces our modified 
equation (2.13). (This trick was also used in [13] in the special case 9 — 0.) 


Remarks to mathematicians. The above argument contains some sleight of hand. 
In particular, the domain 0 was replaced by the halfspace s > 0 or, equivalently, a > 0. 
This change removes any ordinary boundary layer which may be present at the outflow. 
In addition, boundary data for (2.13) are to be applied at the rotated boundary s = 0. 
We expect these distortions to introduce discrepancies only near the point of discontinuity 
(x,j/) = (0,j/o)- In particular, it appears from our computations that there needs to be 
a slight shift of coordinates. See the comments concerning Fig. 1 below. We have also 
not shown that the solution of the modified equation (2.13) bears any resemblance to the 
solution of the difference scheme (2.8). Such a proof would probably proceed as in [13] with 
the replacement of 0 by the halfspace x > 0, followed by a Fourier transformation of (2.8) in 
the y-direction. The modified equation (2.13) shows that the canonical form of the integral 
representation of v is 


v(<t,t) 


+ 


i: 


/(u>) f. 3 

exp \ UI 3 UJ 

*■ 


a 2 w + 


| dui 


( 2 . 21 ) 


with / and the cij dependent on <r and r. Furthermore, the integral (2.16) indicates that 
/ « 1, a 3 « jSc, a 2 r^> (T , and dj ~ r when cr 1 and |r| 1. The integral (2.16) derived 

from the modified equation (2.14) is merely a nonuniform asymptotic approximation which 
is valid when |r[ and when (3 is bounded away from zero. We see from the form 

of (2.21) that a uniform asymptotic estimate would require investigation of the interaction of 
two saddle points and a pole. For the case when sin(40) = 0 the situation is simpler because 
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<*3 — 0 and there is only one saddle point and a pole. Uniform asymptotics for 9 = 0 are 
presented in Hedstrom and Osterheld [13]. 

A computational example. In our computations to illustrate these oscillations we 
located the point of discontinuity at y 0 = 0.25, we chose coefficients 


F cos 0 = 2, Fsin0 = l, e = 0.002, 


and we used a mesh size of h = 0.02. This gives a cell Reynolds number of moderate size 
Rh = 10\/5, and with L = 1 it gives 7 = 0.02. The scaling (2.10) is therefore \J~Ltjv & 
0.0946, and the value of /? in (2.15) is /3 fa 0.598. The cross section at * = 0.8 is shown in 
Fig. 1, where the solution to (2.8) is shown as a solid curve and the Airy integral (2.16) is 
given as dashes. We must admit that in order to obtain such a good match of the curves, 
we had to shift the jump for the Airy integral from y 0 to y 0 + h. This could be because the 
Airy integral applies to the rotated coordinate system (s,t) given by (2.9). It shoul<J also be 
noted that there is a phase difference between the two curves in the oscillatory region. This 
is a well-known deficiency of modified equations, and it results from the nonuniformity of 
the asymptotic approximation. At the point (x,y) = (0.8, 0.6) near the overshoot the value 
of the parameter a given by (2.17) is a « 1.294. We have oscillations because a < 2. 

The numerical method we used to solve (2.8) is a combination of ideas from Elman 
and Golub [9] and from Chin and Manteuffel [6]. As in Elman and Golub, we introduce 
a red-black ordering on the grid points and do a cyclic reduction to obtain a nine-point 
scheme on the black grid points. This reduction produces a matrix much better conditioned 
for iterative methods. The iterative method used by Elman and Golub is point Jacobi, 
mostly because they impose no constraints on the direction of flow. In our example the 
flow is one- directional, so we follow Chin and Manteuffel in using line Gauss-Seidel with 
lines transversal to the flow, starting at the inflow boundary and marching downstream. We 
find that this scheme converges very rapidly, with the greatest speeds at high cell Reynolds 
numbers. (Perhaps, we should reiterate that the point of this section is to show that rapid 
solution of the matrix equation should not be the primary objective — its solution is a poor 
approximation to the solution of the differential equation when the parameter /? in (2.15) is 
large.) 


Let us remark that we have also solved (2.8) in a version with a discrete approximation 
to Neumann outflow boundary conditions 


d x u = 0 for x = 1, 
d y u = 0 for y = 1. 


( 2 . 22 ) 


We found this boundary condition to be satisfactory only for small cell Reynolds number, 
Rh < 5. Otherwise, there are additional small oscillations with period 2 h induced by the 
mismatch at the outflow boundary x = 1. 
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3 Curvilinear coordinates. 


In this section we permit the coefficients a and b in (1.1) to depend on the position (jc,y), 
and we present a numerical algorithm for generating an orthogonal coordinate system (a 
chart) aligned with the given vector field (a, b). Our coordinate system is derived from the 
characteristic curves. We remark that a somewhat different coordinate transformation based 
on characteristics was given by Chin et aJ. [5]. 

We again assume that the vector field (a, b) has no stagnation point, so that |a| + |6| is 
bounded away from zero for all (x,j/) in ft. For purposes of constructing the mapping, it 
is convenient to do an initial scaling so that a 2 + b 2 — 1. One of our goals is to set up a 
mapping (s,f) i— > (x,y) such that s follows the flow in the sense that there exists a positive 
function cf> for which 

d, = <j){a d x + bd y ). (3.1) 

Because the vector (—b, a) is orthogonal to (a, 6), the orthogonality requirement (our second 
goal) amounts to the condition 

d t = if>{-bd x + ad y ) (3.2) 

for some positive function In a moment we show that the scale factors <f> and i}> are not 
arbitrary. 


In part, the construction of such a mapping is easy, because it is easy to integrate (3.1). 
All that is needed is to pick a convenient starting point (* 0 , 3 / 0 ) and to integrate the system 


dx 

ds 

dy 

ds 


a(f>, 

bd>, 


x — xq at s = 0, 
y = y 0 at s = 0. 


(3.3) 


This gives a curvilinear coordinate line in ft corresponding to a constant value of f. The 
image of a line s = constant may be obtained similarly by integrating 


dx 

It 


= -H, 



x = xq at t = 0, 
y = y 0 at t = 0. 


(3.4) 


We still must ensure global consistency as follows. Let us traverse the edges 'of the 
curvilinear rectangle $o < s < s\, to < t < t\, and we assume that this rectangle is contained 
in ft. Denote the image of (so,^o) as the vertex A. Suppose further that we integrate (3.3) 
from s 0 to Sj, arriving at the vertex C. We then integrate (3.4) from t 0 to t\ and arrive 
at the vertex B opposite A. Let us now reverse the order by first integrating (3.4) from to 
to ti to arrive at the vertex D and then integrate (3.3) from so to sj. Can we be certain 
that we again arrive at the vertex B? It happens that this global consistency question has 
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been answered [15], and that what is required is the vanishing of the Lie bracket = 

d,d t - d t d,. 


It is easy to see by a short computation that the vanishing of the Lie bracket [$,,(?*] is 
equivalent to the system of partial differential equations 


= d t {b<j>), 

d,(bip) = -d t (a<f>). 


(3.5) 


Upon differentiating the products in (3.5) and solving for d,ti> and d t 4>, we find that a 
necessary and sufficient condition for consistency is that 


(a 2 + b 2 ) d,ij> = <f>(a d t b — bd t a ) — ij>(ad t a + bd,b), 
(a 2 -f b 2 )dt<t> — —<f>(ad t a + bd t b) + ^(bd,a — ad,b). 


(3.6) 


Note that if (a, b) has been scaled so that a 2 + b 2 = 1, then (3.6) takes the simpler form 

(3.7) 


d,Tf> = 4>(ad t b — bd t a), 
d t (f> = ij>{bd t a — ad,b). 


We recognize the system (3.7) as the telegraphers’ equation, written in terms of Lie 
derivatives along the characteristic curves. Therefore, all that is needed for its solution is to 
prescribe values <f> = 1 at t = 0 and = 1 at s = 0 and to march in the s and f-directions 
concurrently. 

It should be emphasized that theoretical questions remain for this grid-generation scheme. 
In particular, there is no guarantee that the solutions </> and ^ will be positive at all points 
in 0. This is important in that the Jacobian of the transformation (3.3-4) is given by 
J = ( a 2 + 6 2 )c^>. We required at the outset that a 2 -f b 2 be bounded away from zero. Thus, if 
we are to maintain a nonzero Jacobian, we must take special measures whenever it happens 
that <f> < 0 or $ < 0. One possibility is to back up and put a boundary on this local chart. 
We could then initialize a new chart and continue. 


4 Domain decomposition for an internal layer. 

In this section we present a computational example which uses domain decomposition to 
resolve an internal layer. At this point we have not yet implemented the algorithm described 
here, but the final report will have computations. In our algorithm we first identify the 
internal and boundary layers, and we then set up a domain decomposition to segregate 
them. The domain decomposition is carried out with overlapping grids using the tools of 
Chesshire and Henshaw [2], We have added the modification that in some subdomains we 
use the grid- generation algorithm of Section 3. 
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As our domain 0 we use the square 0<x<l,0<y<l, and on fi we consider the 
convection-diffusion equation 

(1 -f x)d x u + (1 - y)d y u - eAu. (4.1) 

As boundary conditions for (4.1) we prescribe u = 0 on the bottom of (y = 0), u — 1 on 
the left-hand edge (x = 0), u = 1 on the top ( y = 1), and u = — 1 on the right-hand edge 
(x = l). 

Note that in (4.1) we have chosen coefficients so that there is no turning point in 12. That 
is, we have |1 + x| -f |1 — y\ bounded away from zero in fl. Note also that by the discussion in 
Section 1 the inflow boundary T/ consists of the bottom y — 0 and the left-hand side x = 0 
of the square ft. Furthermore, the top of the square y = 1 is a tangential boundary Ty, and 
the right-hand edge x = 1 is an outflow To- The reduced equation is 

(1 + x)d x TJ + (1 - y)d v U = 0, (4.2) 

and its boundary conditions are imposed on the inflow boundary T/. It so happens that 
we can write down a formula for the solution U of (4.2), although this is not necessary for 
our domain-decomposition algorithm. The characteristic curves for (4.2) are the hyperbolas 
(x — l)(t/ + 1) = const. Thus, the solution of the reduced equation (4.2) is 

tj _ f 1 if V > x/(x + 1), 

(0 if y < x/(x + 1). 

This gives us an internal layer along the hyperbola y = x/(x -f 1) and exponential boundary 
layers at the outflow boundary x — 1. It happens that we imposed boundary data along 
the tangential boundary Tr such that no boundary layer resides there. If there had been a 
boundary layer along IV, we could have modified the domain decomposition described below 
so as to include its effects. 

As the problem is stated, we need the following subdomains: (1) a square B of diameter 
0(e) at the origin to cover the birth of the internal layer, (2) an internal-layer region 

1 = {(x,y): | y - x/(x + 1)| < Gy/tx} 

with 0(e) < x < 1 — 0(e), (3) three outflow boundary layers O, one above the internal layer, 
one below it, and one interacting with it, (4) an outer region Ti above the internal layer on 
which uwl, and (5) an outer region 7i below the internal layer on which u « 0. 

In the two outer regions T~C we use a coordinate system derived from the characteristics, as 
described in Section 3. In the internal layer J we use a parabolic coordinate system imposed 
on the characteristics. (More precise details will be given in the final report.) Finally, in the 
birth B and boundary-layer regions O we use the methods given in the papers by Hedstrom 
and Howes [11] and [12]. The iterations are performed in the order: (1) the outer regions 
(2) the birth region B, (3) the internal layer T, (4) the outflow boundary layers O. The 
iterative schemes in the subdomains are as in [11] and [12]. 
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